maxpop = max(c(out_model_punkty_adresowe$TOT, out_model_punkty_adresowe$EST_POP))
p1 <- ggplot(data = out_model_punkty_adresowe) +
geom_sf(aes(fill = TOT)) +
scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop)) +
labs(title = "Ryc.1 Mapa oryginalnych wartości liczby ludności") +
theme_bw() + theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
p2 <- ggplot(data = out_model_punkty_adresowe) +
geom_sf(aes(fill = EST_POP)) +
scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop)) +
labs(title = "Ryc.2 Mapa estymowanych wartości liczby ludności") +
theme_bw() + theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
print(p1)
print(p2)
p3 <- ggplot(data = out_model_punkty_adresowe) +
geom_sf(aes(fill = RES)) +
scale_fill_gradient2(name = "Błąd estymacji", low = "blue", high = "red", limits = c(min(out_model_punkty_adresowe$RES), max(out_model_punkty_adresowe$RES))) +
labs(title = "Ryc.3 Mapa różnic między wartością obserwowaną liczby ludności a jej estymacją") +
theme_bw() + theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
print(p3)
out_model_punkty_adresowe$RES_reclass = case_when(
out_model_punkty_adresowe$RES < 0 ~ "bledy ujemne",
out_model_punkty_adresowe$RES == 0 ~ "brak bledu",
out_model_punkty_adresowe$RES > 0 ~ "bledy dodatnie",
)
p4 <- ggplot(data = out_model_punkty_adresowe) +
geom_sf(aes(fill = RES_reclass)) +
scale_fill_manual(name = "Typ błędu",
values = c("bledy dodatnie" = "red", "brak bledu" = "white", "bledy ujemne" = "blue")) +
labs(title = "Ryc.4 Mapa zreklasyfikowanych błędów modelu") +
theme_bw() + theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank())
print(p4)
meanerror = mean(out_model_punkty_adresowe$RES)
rmse = sqrt(mean((out_model_punkty_adresowe$RES)^2))
print(paste("średni błąd estymacji modelu wynosi", round(meanerror,2), "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))
## [1] "średni błąd estymacji modelu wynosi 3.67 natomiast pierwiastek średniego błędu kwadratowego wynosi 30.44"
Maciejo,
YOU KNOW WHAT TO DO
# wykres rozrzutu
ggplot(pop1km, aes(x=TOT, y=N_BUDYNKI)) + geom_point() +
coord_fixed(ratio = 1) +
xlab("Liczba ludności") +
ylab("Liczba mieszkań") +
geom_smooth() +
geom_abline(color = 'grey') +
theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
# model liniowy
model_liniowy <- lm(TOT ~ N_BUDYNKI, data = pop1km)
plot(model_liniowy)
W jakim stopniu liczba ludności (TOT) jest wyjaśniana przez liczbę mieszkań?
resid_interact(model_liniowy)
## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## ℹ The deprecated feature was likely used in the ggplot2 package.
## Please report the issue at <https://github.com/tidyverse/ggplot2/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The following aesthetics were dropped during statistical transformation: label
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
summary(model_liniowy$residuals)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -98.6520 -2.5883 -0.7094 0.0000 -0.5883 174.7735
mean(model_liniowy$residuals)
## [1] 3.738246e-16
Wartości odstające można zidentyfikować wykorzystując wykresy diagnostyczne, a w szczególności wykres Residuals vs. Leverage.
autoplot(model_liniowy) + theme_classic()
#dopasowanie modelu
smr_model_lm <-summary(model_liniowy)
c( R2 = smr_model_lm$r.squared, R2_adj = smr_model_lm$adj.r.squared)
## R2 R2_adj
## 0.7019419 0.6990481
#~70% zmienności zmiennej TOT jest wyjaśniona przez zmienność zmiennej N_BUDYNKI.
# Określenie dowolnego przedziału ufności dla współczynników modelu
confint(model_liniowy, level = 0.99)
## 0.5 % 99.5 %
## (Intercept) -7.022728 8.441527
## N_BUDYNKI 3.275627 4.603257
outlierTest(model_liniowy)
## rstudent unadjusted p-value Bonferroni p
## 21 9.742598 3.0248e-16 3.1760e-14
## 25 6.813338 6.7995e-10 7.1395e-08
## 34 -5.176442 1.1389e-06 1.1958e-04
ggplot(pop1km, aes(x=TOT, y=N_BUDYNKI)) +
geom_point() +
geom_point(data = pop1km[c(21, 25, 34),], aes(x=TOT, y=N_BUDYNKI), color = "red", size = 4) +
coord_fixed(ratio = 1) +
xlab("Liczba ludności") +
ylab("Liczba mieszkań") +
geom_smooth() +
geom_abline(color = 'grey') +
theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
par(mfrow = c(2, 2))
plot(model_liniowy)
meanerror = mean(model_liniowy$residuals)
rmse = sqrt(mean((model_liniowy$residuals)^2))
print(paste("średni błąd estymacji modelu wynosi", meanerror, "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))
## [1] "średni błąd estymacji modelu wynosi 3.73824648485021e-16 natomiast pierwiastek średniego błędu kwadratowego wynosi 26.25"
PYTANIA Czego dowiadujemy się z wykresów diagnostycznych? Czy model jest dobrze dopasowany? Czy spełnione są założenia regresji liniowej?
# Usunięcie wierszy 21, 25 i 34
pop1km_n <- pop1km[-c(21, 25, 34), ]
# wykres rozrzutu
ggplot(pop1km_n, aes(x=TOT, y=N_BUDYNKI)) + geom_point() +
coord_fixed(ratio = 1) +
xlab("Liczba ludności") +
ylab("Liczba mieszkań") +
geom_smooth() +
geom_abline(color = 'grey') +
theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
# model liniowy
model_liniowy_n <- lm(TOT ~ N_BUDYNKI, data = pop1km_n)
autoplot(model_liniowy_n) + theme_classic()
#dopasowanie modelu
smr_model_lm_n <-summary(model_liniowy_n)
c( R2 = smr_model_lm_n$r.squared, R2_adj = smr_model_lm_n$adj.r.squared)
## R2 R2_adj
## 0.8914342 0.8903486
#~89% zmienności zmiennej TOT jest wyjaśniona przez zmienność zmiennej N_BUDYNKI.
# Określenie dowolnego przedziału ufności dla współczynników modelu
confint(model_liniowy_n, level = 0.99)
## 0.5 % 99.5 %
## (Intercept) -2.413594 3.198192
## N_BUDYNKI 3.198037 3.843294
outlierTest(model_liniowy_n)
## rstudent unadjusted p-value Bonferroni p
## 65 -6.324628 7.3790e-09 7.5266e-07
## 33 4.742867 7.0854e-06 7.2271e-04
ggplot(pop1km_n, aes(x=TOT, y=N_BUDYNKI)) +
geom_point() +
geom_point(data = pop1km_n[c(65, 33),], aes(x=TOT, y=N_BUDYNKI), color = "red", size = 4) +
coord_fixed(ratio = 1) +
xlab("Liczba ludności") +
ylab("Liczba mieszkań") +
geom_smooth() +
geom_abline(color = 'grey') +
theme_bw()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
#WIZUALIZACJA
par(mfrow = c(2, 4))
plot(model_liniowy)
plot(model_liniowy_n)
meanerror = mean(model_liniowy_n$residuals)
rmse = sqrt(mean((model_liniowy_n$residuals)^2))
print(paste("średni błąd estymacji modelu wynosi", meanerror, "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))
## [1] "średni błąd estymacji modelu wynosi 2.30805383656037e-16 natomiast pierwiastek średniego błędu kwadratowego wynosi 9.14"
ggplot() +
geom_sf(data = pop1km, aes(fill = TOT)) +
scale_fill_viridis_c()
labs(fill = "Liczba ludności", title = "Mapa liczby ludności") +
theme_minimal()
## NULL
# Przewidywanie EST_POP na podstawie modelu liniowego
pop1km$EST_POP <- predict(model_liniowy, newdata = pop1km)
# Tworzenie mapy z estymowaną liczbą ludności (EST_POP)
ggplot() +
geom_sf(data = pop1km, aes(fill = EST_POP)) +
scale_fill_viridis_c() + # Możesz zmienić paletę kolorów wg potrzeb
labs(fill = "Estymowana liczba ludności", title = "Mapa estymowanej liczby ludności") +
theme_minimal()
# Obliczenie reszt
pop1km$Reszty <- pop1km$TOT - pop1km$EST_POP
# Tworzenie mapy reszt
ggplot() +
geom_sf(data = pop1km, aes(fill = Reszty)) +
scale_fill_viridis_c() + # Możesz zmienić paletę kolorów wg potrzeb
labs(fill = "Reszty", title = "Mapa reszt (TOT - EST_POP)") +
theme_minimal()
# Przeklasyfikowanie reszt do trzech klas
pop1km$Klasy_reszt <- cut(pop1km$Reszty, breaks = c(-Inf, 0, Inf), labels = c("Przeszacowane", "Niedoszacowane"))
# Utworzenie etykiety dla braku błędu
pop1km$Klasy_reszt[is.na(pop1km$Klasy_reszt)] <- "Brak błędu"
## Warning in `[<-.factor`(`*tmp*`, is.na(pop1km$Klasy_reszt), value =
## structure(c(1L, : niepoprawny poziom czynnika, wygenerowano wartość NA
# Tworzenie mapy przeklasyfikowanych reszt
ggplot() +
geom_sf(data = pop1km, aes(fill = Klasy_reszt)) +
scale_fill_manual(values = c("Przeszacowane" = "red", "Brak błędu" = "grey", "Niedoszacowane" = "blue")) + # Możesz zmienić kolory wg potrzeb
labs(fill = "Klasy reszt", title = "Mapa przeklasyfikowanych reszt") +
theme_minimal()
# Tworzenie interaktywnej mapy rozmieszczenia ludności w siatce 100m
pop_100m_plot <- ggplot() +
geom_sf(data = pop_100m, aes(fill = TOT)) +
scale_fill_viridis_c() +
labs(fill = "Liczba ludności", title = "Mapa rozmieszczenia ludności w siatce 100m") +
theme_minimal()
# Konwersja ggplot na obiekt plotly
pop_100m_plotly <- ggplotly(pop_100m_plot)
# Wyświetlenie interaktywnej mapy
pop_100m_plotly